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Introduction 

Composite systems constitute a large class of naturally occurring or artificially 
synthesized disordered systems [p]]. The systems are microscopically inhomogeneous 
and disordered but look homogeneous on the macroscopic scale. From the tunnelling 
electron micrographs (TEM) of such a composite material it can be seen that the 
typical dimension (£) of metallic islands embedded in the insulating matrix are much 
greater than the atomic size (a) but obviously much smaller than the macroscopic 
scale length (L): a <C £ <C L. The effective conductivity of such a system depends 
upon the conductivities of the individual phases. For a low volume fraction (p) of the 
conducting phase, the system as a whole behaves like an insulator since the conducting 
regions do not form a continuous path through the sample. As p is increased, the 
conducting regions will in general tend to grow and eventually at a critical volume 
fraction (p c , called the percolation threshold) the conducting phase percolates through 
the sample. This may be considered as a classical insulator-to-metal transition or 
more popularly as a percolation transition. For all p > p c , the system is metallic, and 
if the conducting phase is Ohmic, so is the whole macroscopic system. Clearly this 
class of systems may be well described by the geometrical percolation theory. 

Now if an external voltage is applied across such composite systems (examples 
include dispersed metallic systems, carbon-black-polymer composites, sulphonated 
(doped) polyaniline networks etc., which are usually highly structured and give rise 
to some sort of universal behaviour.) a wide variety of interesting features asso- 
ciated with a nonlinear response emerge. Usually these composites exhibit an un- 
usually low percolation threshold. Qualitatively identical nonlinear I — V (as well as 
dl/dV (= Gjagainst V) response have been reported [Q, [| both below and above 
the threshold for many of the composites although the nonlinearity exponent is found 
to be grossly different in the two regimes. Power-law growth of excess conductance 
for small V is another general feature of the class of composite systems where non- 
integer power-law has been observed. This in turn implies a power-law in the I — V 
relationship for small applied voltage (V). The G — V curves are seen to saturate 
for an appropriately high enough voltage below the Joule-heating regime. The typ- 
ical curve then looks like a nonlinear sigmoidal type function interpolating two lin- 
ear regimes. Recent experiments on carbon-wax systems [fj] as well as many earlier 
ones on disordered/ amorphous systems [|J, find a non-integer power-law behaviour 
and a saturation in the DC-response as mentioned above. Composite systems show 
very interesting temperature-dependent conduction properties particularly in the low 
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temperature regime where the conduction is mainly due to phonon- assisted hopping 
(Mott's variable range hopping (VRH)). Some recent experiments indicate an effect 
of dilution on the relevant temperature-exponent for fitting the low-temperature data 
with VRH or its other variations. We address many of the above mentioned features 
through our study. 

Modelling Nonlinear Transport 

The framework of our study is based on percolation theory. The ultra-low perco- 
lation threshold and the fact that many of these nonlinear systems carry current even 
below p c indicates strongly that tunnelling through disconnected (dispersed) metallic 
regions must give some virtually connected percolating clusters. From the nonlin- 
ear I — V characteristics (e.g., see the experiment by Chen and Johnson ||) it is 
observed that the response (DC) behaviour is reversible with respect to the applied 
field. Also the temperature-dependent resistance with a minimum at some charac- 
teristic temperature and the Mott variable range hopping (VRH) type behaviour at 
very low temperatures give further credence to tunnelling assisted percolation. In 
practice, the tunnelling conductance should fall off exponentially with distance and 
hence the tunnelling should have an upper cut-off length scale. So for simplicity and 
to capture the basic physics, we construct a bond (lattice) percolation model for this 
problem, such that tunnelling may take place only between two nearest neighbour 
ohmic bonds and no further. For a further simplification, we assume the nonlinear 
response of each tunnelling bond (or resistor) to be piecewise linear. We assume that 
all the tunnelling bonds have an identical voltage threshold (v g ) below which they are 
perfect insulators and above which they behave as ohmic conductors. Clearly this is 
the source of nonlinearity in the model. Made of both random resistive and tunnelling 
elements, this network will be called a random resistor cum tunnelling-bond network 
(RRTN) Q. 

The percolation statistics || of the model network is examined in the saturation 
limit, i.e., when all the tunnelling bonds can overcome their threshold. We estimate 
the new percolation (p ct ) threshold and address the question of universality class. We 
undertake small-cell renormalization, Monte Carlo simulation and finite size scaling 
analysis to estimate p ct and some of the idependent critical exponents around it. The 
simulation results are obtained for lattices in 2D for convenience. Lattice sizes L = 20 
to 300 are considered. The p ct is found to be 0.181 ± 0.001 and the value of correlation 
length exponent v is obtained to be = 1.35 ± 0.06. The fractal dimension (D) for 
the spanning cluster at the threshold which is found to be D = 1.87, very close to 
91/48, the fractal dimension for 2D random bond percolation. We also calculate the 
conductivity exponent, t, in the upper linear regime where all the tunnelling resistors 

1 In this respect we comment that a dynamic random resistor (DRRN) model proposed by Gefen 
et at, is different from our RRTN model in the sense that they allowed any insulating bond 
at any position in the lattice to break and turn metallic, whereas in our case such breakings can 
occur only at some correlated bond positions. Moreover, with the addition of these bridge bonds 
(anywhere), the new percolation threshold may not be properly defined like that of ours. 
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are considered to be behaving as the other ohmic resistors. We obtain tjv = 0.90, 
where the value of this ratio for the usual percolation problem is = 0.97. The value 
of the critical exponents, as obtained above, indicate that this correlated model for 
percolation belongs to the same universality class as that of its uncorrelated version 
(in the absence of tunnelling bonds). 

An effective medium approximation (EMA) |J has been used to calculate 
the percolation thereshold for our model system and the conductivity behaviour in 
the saturation limit. The probability of a bond to be ohmic, tunnelling or purely 
insulating according to the considerations of our model is: P hm = P, Ptun = (p 3 + 
3p 2 q + 3pq 2 ) 2 q, Pi ns = [1 — (p 3 + 3p 2 q + 3pq 2 ) 2 ]q, where q—l—p. If the conductances 
of the three types of bonds are denoted by g h m , Qtun and gi ns , the EMA equation for 
this general situation can be written as 

Pohmifle Qohm) Ptun{G e Qtun) PinsiGe Qins) q 

[9ohm + (d-l)G e ] [g tU n + (d-l)G e ] \g~ + (d-l)G e }~ 

Solving the above equation for 2D (d = 2) and 3D (d = 3) we obtain p ct = 1/4 and 1/8 
respectively. It may be noted that the value of p ct in 2D is close to that obtained by 
numerical simulation. We examine the behaviour of the effective linear conductance 
(G e ) for the macroscopic model composite system in the saturation limit, given some 
specific values or forms of the resistive elements. These then are compared with the 
results obtained with the numerical simulation. The agreement is fairly good when 
one is away from the threshold, p ct . 

Next we study the dielectric breakdown phenomenon in our model as the on- 
set of nonlinear conduction against applied field for p < p c . Below the percolation 
threshold (p c ) there exists a number of metallic clusters, isolated from each other, 
but closely spaced. As new conducting paths are created when the local electric 
field across tunnelling bonds increases above v g , the conductivity of the whole system 
jumps from a zero to a non-zero value (for p < p c ) as the external applied field crosses 
[0] the dielectric breakdown field (E B = V B /L). Note however that below p ct there 
is no sample-spanning cluster of combined ohmic and tunnelling bonds, and hence 
there is no breakdown at any finite electric field according to the criterion set for 
our model. The interest would be to estimate the breakdown exponent t B , where 
E P B ~ (p c — p) tB . To remove finite size effects, we work with the asymptotic break- 
down field E B (L = oo). From the least-square fit of the data for the above we find 
that the breakdown exponent ts — 1-42 for our RRTN model. It seems that the 
above exponent is not very different from that of the usual breakdown exponent 
tg = v =1.33 as discussed above. But it is not unlikely either that we have a different 
result in our hands. If different, it could be because of the nature of the electric field 
in increasing the effective volume fraction of the conductors. 

We present the nonlinear DC-response namely, the current- voltage (I — V) 
and the conductance-voltage G — V charactersistics in our model system. Our com- 
puter simulation involves solving Kirchhoff 's law of current at the nodes of the RRTN 
network in 2D with the linear and nonlinear (assumed piecewise linear) resistors and 



3 



20 
18 
16 
— 14 
r 12 
10 



3 
O 



8 
6 
4 
2 




"i 1 1 r 

L=20 



"i r 



i r 



p=0.9 Q □ 



□ 



□ -I 



@ p£> <> ' ^ — i 1 1 1 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



□ 



+ 



p=0.5 + + + 
+ p=0.3 



+ 



+ 



+ 



o o 



o o 



o o 



o o 



o o * 



2 4 6 8 10 12 14 16 
Applied Voltage (V) 



18 20 



Figure 1: Current (J) against Voltage (V) curves for different volume fractions (p) of 
the conducting components. 



the standard Gauss-Seidel relaxation technique. We obtain current (I) and therefrom 
the differential conductance (G = dl/dV) for the whole network at a given volume 
fraction p of the ohmic bonds. Simulation results for nonlinear I — V curves for a 
square network of size L = 20 were plotted in fig. 1 for p = 0.3, 0.5 and 0.9. Averages 
over 50 configurations are done in each case. One may note that the nonlinearity in 
the response exists for all p both below and above p c . 

The differential conductance G (= dl/dV) of the network is obtained directly 
from the I — V curves. A typical such G — V curve is shown in fig. 2 for L = 20 and 
p = 0.8. To understand the conductance behaviour for the entire network we adopt 
a pedagogical approach where we analyse the elementary prototype circuits with 
nonlinear resistors. The conductance (G) of these elementary units grows nonlinearly 
with the applied voltage V and gives us an idea of what type of functions may be used 
to fit the G — V data for the much more complex macroscopic system. After sifting 
through various such functional froms, we find that the simulation data obtained 
through our model system in 2D were best fitted with: 

G = G + G d [l - exp(-\AW)]\ (2) 

where Gd = Gf — G and AV = V — V g , where V g is the macroscopic threshold 
voltage and is the same as Vb above. Go is the conductance in the limit AV — > 0. 
Experimentally Gf may be obtained by applying a large enough voltage (V s ) such 
that Joule heating remains unimportant. In our computer simulation on finite sized 
systems, we find V s to be many orders of magnitude larger than V and Gf is the 
conductance when all the tunnelling bonds take part in the conduction. 

For a meaningful comparison of all the G — V data with different Go, Gf, V g , etc., 
we scale the conductance G as G = (G— Go) / Gd and the voltage V as V = (V— V g ) /V g . 
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Figure 2: A typical curve showing the behaviour of differential conductance G against 
V. 



In fact, we tried to scale the G — V data for a set of p in the range 0.48 < p < 0.52 
(i.e., both below and above p c ), and we found that all the data do reasonably collapse. 
This suggets the following general form for the functional behaviour close to p c ; 



G = f(V), 



(3) 



where f(x) is a function such that /(0) = 0, and /(oo) = 1. Here we point out 
that the threshold (or the breakdown) voltage V g is the only relevant voltage scale 
that enters into the scaling function. The other voltage scale V s is seen to have no 
role in the above scaling eqn. (^). Expanding eqn. (|2|) near the onset of nonlinearity 
(AV — > 0), the excess conductance AG = G — G varies with the voltage difference 
(AV) as a power-law: 

AG ~ AV^ = AV 5 , (4) 

and the nonlinearity exponent 5 = /ry. For p close to p c (0.48 < p < 0.52) we find 
that 5 = 1.0. Thus the nonlinearity exponent for the I — V curve is a = 5 + 1 = 2.0. 
Experiments in 2D arrays of normal metal islands connected by small tunnel junctions 
by Rimberg et al. [0] found a = 1.80 ±0.16; suggesting a good support for our model. 

Our analysis of results for widely different volume fractions indicate that the 
nonlinearity exponent 6 increases significantly as we go sufficiently away (both below 
and above) from the percolation threshold. The scaled data for all the curves now 
do not fall on top of each other indicating that all of them can not be described by 
the same fitting function f(x) or by the same fitting parameters fi and 7. Hence the 
possible power-law in the regime (AV — ► 0) for the onset of nonlinearity for these 
curves of different p are not all the same. 

Now the difference between two limiting conductances, Gd = G/ — Go may be 
taken as a measure of overall nonlinearity, whereas the nonlinearity exponent (S or 
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a) gives a measure of the initial nonlinearity near the threshold. Gd as a function of p 
shows a peak at around p = p c . So we find that the overall nonlinearity is maximum 
near the geometrical percolation threshold. Next we looked at how Gd is related to 
the initial conductance Go i n the interval p ct < p < 1 in the limit L — > oo. The 
relationship is linear which actually means that Gf is also linearly dependent on Go. 
This in turn implies an identical p-dependence for the two saturation conductances 
G and Gf around their respective thresholds (p c and p ct ), consistent with the fact 
that the system has the same conductivity exponent in both the zero and the infinite 
voltage limits (i.e., G j ~ (p -p c ,ctY)- 

The AC-response of the model system also turns out to be very interesting. 
In this case the tunnelling bonds in RRTN are assumed to behave as capacitors. 
The AC-conductance is now expected to behave nonlinearly between two saturation 
regions of u — > + and u> — ■> oo as in the case of DC-response discussed above. We 
first give here the EMA where each tunnelling bond has the conductance g tun = iwc, 
where i = a/— 1, c is the capacitance of the tunnelling bonds and uj is the circular 
frequency of the applied sinusoidal voltage with unit amplitude. Here we take c = 1 
for convenience, thereby setting the frequency scale. So for a square lattice (d = 2) if 
we take g hm = 1 in eqn. (|l|), the real part of G e (u) can be shown to be 

ReG e (co) = (2Pofe ™ " 1] + i(X 2 + Y^cos 6 -, (5) 

where X = [2P ohm - l) 2 - u\2P tun - if and Y = 2u[(2P ohm - l){2P tw - 1) - 
2(2P ins — 1)] and 9 = tan _1 (y/X). It may be checked from the above eqn. (|5|) that 
at p = p c (= 1/2 in 2D) and in the limit u; — ^ the real part of the complex effective 
conductance behaves as ReG e (uj) ~ u 0,5 . This is also true for 3D. 

Next we look at the simulation results for the AC-response. It has been observed 
that for frequencies u < u , one gets some generic linear or quadratic dependences on 
to which may be easily understood. But, for frequencies to > u , we expect percolative 
effects to gain control and G rma (u) to follow an equation similar in form to that used 
for the DC-conductance: 

G rms (uj) = G rms (uj ) + G d (u)[l - exp(-A[w - ^oHF- (6) 

For many practical situations, the intermediate frequency range (between uq and 
the upper saturation) is of prime interest. In this case, fitting the average graphs, 
we find that S'(= fij) has a minimum value of about 0.7 near p c , and increases on 
both sides of it. In other words, the AC nonlinearity exponent 5' (away from p c ) is 
also p-dependent. Notwithstanding this fact, experiments |, | on a wide variety of 
disordered systems observe 5' = 0.7 which matches closely with our result. 

The behaviour of phase-angle (sometimes called the loss-angle) of the complex 
conductance G e (u) with respect to frequency (lu) is of practical interest. The phase- 
angle (0) is defined through tan0 = ImG e (u)/ReG e (u). The shift of the peak value 
of it with the dilution (p) is worth noting. The agreement between the simulation 
result and that by EMA is reasonably good. The variation of the phase-angle (0) 
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with frequency (u) has been observed for a range of values 0.3 < p < 0.7. The peak 
value cj) m = 0.7 (radian) is obtained for p around p c which is close to the universal 
phase-angle value of 7r/4 obtained in the simple RC model in 2D at p c predicted by 
Clerc et. al. [|J. We looked at the phase-angle versus frequency for p = p c , and we 
find therefrom that 4> m = 0.7 in this model too. 

Composite systems have very interesting temperature dependent conduction 
properties |9j particularly in the low-temperature regime. Some recent experiments 
on them show the analysis of their low-temperature data which seem to be confusing 
and contradicting each other. The controversy, as briefly described below, is still 
on and the complete physics is yet to be understood. The usual attempt is to fit 
the low-temperature data for such systems by the well-known Mott variable range 
hopping (VRH) formula or with any of its many generalized forms. In a very recent 
experiment by Reghu et al, in proton-doped polyaniline networks, it was found 
that the exponent in VRH systematically increases from 0.25 to 1 upon decreasing 
the volume fraction p of the conducting component. Here our goal is not to explain 
the recent experimental results exactly. Rather, our modest hope is to demonstrate 
the fact that if one represents the low-temperature data in such systems by the VRH 
or any of its generalized forms, then the the relevant exponent in that can change 
continuously with dilution. The approach is again based on percolation theory where 
we assume the the activated behaviour for the tunnelling bonds and the metallic 
behaviour to the ohmic bonds. The effect of dilution of the temperature dependent 
conductance behaviour can thus be understood at a preliminary level [|T0|. 



Discussion 



In this report we have discussed various aspects of the nonlinear response in the 
disordered binary composite systems in general. We have proposed a very simple 
and minimal model in order to understand the nonlinear electrical response and as- 
sociated physics in composite systems. In many other physical systems, the response 
is negligibly low (or there is no response at all) until and unless the driving force 
exceeds a certain threshold value. So a class of problems exist where sharp thresholds 
to transport occur. The examples in the electrical case is a Zener diode and in the 
fluid permeability problems a Bingham fluid (where there is a critical shear stress r c , 
above which it has a finite viscosity and below which it is so enormously viscous that 
it does not flow). So all these problems may be treated in a similar footing with the 
underlying percolation geometry. 
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